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Abstract 

A  Fuel  Cell  (FC)  is  an  electrochemical  device  which  produces  electric  energy  in  DC.  In  order  to  support  design  control  for  the  electrical  system 
connected  to  it,  it  is  necessary  to  work  out  a  suitable  representation  of  the  fast  dynamics  involved.  Therefore,  in  this  work,  a  mathematical  model, 
based  on  first  principles  and  including  both  dynamical  equations  and  algebraic  relations,  is  described  for  electrochemical  reactions,  with  the  related 
formation  of  potential  differences  and  anion  and  cation  accumulation  phenomena,  in  a  Molten  Carbonate  Fuel  Cell  (MCFC).  The  model  is  formally 
consistent  and  it  has  been  validated  against  experimental  results,  such  as  steady-state  power  and  voltage  versus  current  curves. 

©  2006  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

This  work  focuses  on  electrochemical  reactions  and  on  ionic 
transport  inside  a  Molten  Carbonate  Fuel  Cell  (MCFC),  so  that 
the  consequent  production  of  potential  differences  across  a  cell 
and  the  electric  current  flowing  inside  it  and  through  the  external 
circuit  connected  to  it  can  be  calculated  [1-5].  This  discussion 
can  be  considered  as  a  refinement  of  the  electrical  description, 
based  on  the  Nernst  equation,  presented  in  [6,7],  and  it  is  meant 
to  analyze  relatively  fast  dynamics  (time  scale  around  a  second, 
at  most)  in  the  cell  behaviour.  This  will  help  understand  cell  per¬ 
formance  limits,  especially  when  abrupt  changes  in  the  power 
requested  by  the  user  occur,  as  in  emergency  conditions  (typi¬ 
cally,  either  load  loss  or  sudden  increase  in  power  request).  Of 
course,  this  also  influences  the  design  of  possible  control  sys¬ 
tems,  which,  in  fact,  will  have  to  be  able  to  react  to  such  fast 
and/or  large  “disturbances”  in  a  proper  way,  guaranteeing  power 
quality  and  flexibility  while  keeping  device  operation  safe.  We 
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remark  that,  although,  in  the  literature,  there  are  various  mathe¬ 
matically  exhaustive  models  ([8-10],  e.g.)  for  the  electrochem¬ 
ical  behaviour  of  a  cell,  those  models  hold,  in  general,  for  sta¬ 
tionary  conditions,  i.e.  they  do  not  describe  the  dynamics  of 
the  electrochemical  reactions  at  the  electrodes  and  of  the  ionic 
components  transport. 

2.  The  model 

The  stratified  structure  of  a  single  MCFC  is  reported  in  Fig.  1 . 
Gas  flows  at  the  top  and  at  the  bottom  are  assumed  to  be  in 
orthogonal  directions  with  respect  to  each  other  and  to  the 
cell  thickness:  this  cross  flow,  indeed,  causes  strong  chemical 
coupling  between  corresponding  anode  and  cathode  points  and 
therefore  it  requires  to  study  the  behaviour  of  each  strip  (indi¬ 
cated  with  a  parallelepiped)  along  the  cell  thickness.  The  ap¬ 
proach  followed  here  is  adapted  to  the  three-phase  electrochem¬ 
ical  structure  of  each  strip:  the  gases  (i.e.  the  anodic  and  the 
cathodic  feeding  mixtures)  and  the  liquid  (i.e.  the  electrolyte) 
involved  can  react  only  in  the  presence  of  the  solid  electrode  ma¬ 
terial.  This  last  acts  as  a  catalyst,  offering  empty  sites  to  which 
gas  molecules  can  attach  and  form  activated  complexes,  which, 
in  turn,  react  with  the  electrolyte  ions.  Since  there  are  charged  ac¬ 
tivated  sites,  a  surface  charge  distribution  originates  on  each  cell 
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Fig.  1.  The  considered  MCFC  structure. 

electrode,  so  that  a  potential  difference  arises  and  a  capacitive 
effect  is  highlighted.  That  potential  difference,  together  with  the 
fact  that  carbonate  ions  are  consumed  or  produced  by  reactions, 
is  responsible  of  ionic  movement  inside  the  cell.  The  mathemat¬ 
ical  model  studies  three  separated  zones,  both  at  the  anode  and  at 
the  cathode  (as  summarized  in  Fig.  2,  where  the  spatial  x  coordi¬ 
nate  associated  to  each  of  them  is  shown  too):  the  first,  called  the 
Helmholtz  layer  (compare  [11]),  is  the  one  where  the  electro¬ 
chemical  reactions  take  place  and  it  is  part  of  charged  materials; 
the  second  is  the  electrolyte  zone,  inside  a  pore,  in  which  elec¬ 
troneutrality  cannot  be  assured  (i.e.  it  is  a  space  charge  region), 
since  anions  tend  to  accumulate  near  the  anode  to  decrease  near 
the  cathode;  the  last  one,  which  is  called  here  the  (electrolyte) 
axial  region  and  which  connects  the  anode  and  cathode  cell  sides, 
is  that  in  which  electroneutrality  can  be  assumed. 

2.1.  The  Helmholtz  layers 

The  reaction  mechanism,  described  below  by  formulae  (1)- 
(3)  at  the  anode,  and  (4)— (6)  at  the  cathode,  has  been  split  into 
sub-reaction  steps  which  lead  to  the  formation  of  activa  ted  com¬ 
plexes  (MHC07-  and  MOH-  at  the  anode,  MO-  and  MCO^"_ 
at  the  cathode)  which,  in  turn,  give  rise  to  potential  barriers,  de¬ 
scribed  by  their  peak  values  £pa  and  £ pc .  M  is  a  free  metal  site. 

More  specifically,  the  hypothesized  mechanism  at  the  anode 
is 

1  n 

-H2  +  M^MH,  with  f  1  r_  1  (and  therefore  negligible) 

2  r_i 


anodic  lax  0  cathodic 

radial  x  radial  x 


Fig.  2.  The  adopted  Cartesian  coordinates  and  the  potential  profile  along  a  cell. 


h 

MH  +  CO^“ ^MHCO^“  ^  C02  +  e  +  MOH-  (2) 

r- 2 

rib 

MOH-  +  H2  H20  +  e  +  MH,  (3) 

r-2b 

and  at  the  cathode 

-02  +  M  -M2O2  (4) 

(very  fast,  therefore  negligible:  there  will  always  be  M202  sites) 

1 

-M202  +  e^M0-  (5) 

2  r-  3 

r 4 

MO-  +  C02  +  e^MCO^-  ^  CO^-  +  M.  (6) 

r-4 

r2  and  r_2  are  supposed  to  be  a  fast  and  a  slow  step,  respec¬ 
tively,  rib  a  fast  one  too;  similarly,  r3  and  r_ 4  are  supposed  to 
be  fast  steps,  while  r4  a  slow  one. 

Dynamic  equations  are  based  on  Butler- Volmer-like  reaction 
rates,  and  on  the  Planck-Nernst  formula  for  charged  species 
fluxes.  In  particular,  for  the  anodic  rates  one  can  write 

ri  =  K2C;0(  1  -da-  d.dCt)e-«EP*-E°yRT)F  (7) 

f-2  =  tf-2Cco2a0^e-«£P*-£‘<0+>>/*^  (8) 

rib  =  (9) 

r_2b  =  K_2bCH2 Oaod  -9a-  eact)^ME^0+)-E^RT)F , 

(10) 

where  0a,  @act  and  1  —  9a  —  0act  are  the  surface  fractions  oc¬ 
cupied,  respectively,  by  MOH-,  by  MHCO^-  and  by  the 
H2  molecules  tied  to  the  electrode,  the  Kj’s  are  temperature- 
dependent  coefficients  and  the  concentrations  involved  are  those 
at  the  electrode-electrolyte  interface.  Superscript  (-)  stands  for 
the  carbonate  anion.  According  to  analogous  definitions,  for  the 
cathode  one  can  write 


r3  =  K3Cl£cQ(l  -6C-  0Cact)e-O+feX(^-^(O+))/Sr)F 

(11) 

r_3  =  K-39ce  ^~ERf+)F 

(12) 

r4  =  KhCco2coOcQ  rt  F 

(13) 

£c(0+)-£pC 

r-4  =  K-4Cc0(  1  -Bc-  0Cact)e  rt  f. 

(14) 

Letting  Ta  and  7’act  be  the  number  of  anodic  sites  moles  of  the 
two  possible  kinds  per  unit  area,  the  conservation  equations  of 

the  occupied  sites  are  then 

d  da 

Ta  ,  =  r2-  r-2  -  r2b  +  r_2fc 

dt 

(15) 

d$act 

Tact  dt  =Z  J.d0  (r2  r-2), 

(16) 

848 


S.  Bittanti  et  al.  /  Journal  of  Power  Sources  160  (2006)  846-851 


current  [A] 


> 

<D 

O) 

o 

> 


LU 

I 


o 

o 

HI 


I 

o 

LU 


Fig.  3.  From  top  to  bottom:  power  and  voltage  vs.  current:  simulated  behaviour 
of  a  single  cell  with  700  about  0.8  m2  area;  electrode  potential  differences; 
electrode  surface  fraction  occupied  by  activated  sites,  both  at  the  anode  (i 9a  := 
MOH~,  0act  •=  MHCO^"_)  and  at  the  cathode  ( 6C  :=  MO~,  6C act  :=  MCO^-). 
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where  the  anion  flux  at  the  aO  interface  is,  for  the  Planck-Nernst 
formula, 


•4o — 


_  _  £pa  -  Ea( 0+) 


-v  z  FC. 


a0 


h\/A 


D~ 


d  C~ 

~dR 


(17) 


R  is  the  anodic  radial  x,  /hza  the  Helmholtz  layer  width  and  v~ , 
Z~  and  D~  are  CO3  _  mobility,  number  of  charges  and  diffusion 
coefficient,  respectively. 


Fig.  4.  Gext  from  1500  S  to  1485  S  in  1  s:  from  top  to  bottom:  external  voltage 
and  current,  electrolyte  voltage  and  current,  electrode  potential  differences,  the 
0  variables. 
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time  [sec] 


Fig.  5.  Gext  from  OS  to  1500 S  in  10s:  from  top  to  bottom:  external  voltage 
and  current,  electrolyte  voltage  and  current,  electrode  potential  differences,  the 
6  variables. 


The  anodic  external  current,  ia,  will  be  given  by 

l  Q 

——  =  i~2  -  r_2  +  f2b  ~  r-ib,  (18) 

FSla 

where  5ra  is  the  anodic  reaction  surface.  As  for  the  potential 
difference  Ffl(0+)  —  Ea,  the  classical  Coulomb  relation  supplies 

£fl(0+)-£a  =  -^/hza,  (19) 

^^hza 

with  the  surface  charge  [Cm-2]  in  the  Helmholtz  layer 


°liza  —  Tict ^act ^\] | !('(), —  "h  Fa@aFZy[On  . 


As  before,  for  the  cathode  one  has 

d  0C 

Tc—  =  r3  -  r_ 3  -  r4  +  r_4 
dr 


A6r 


dr 


=  — Z  Tr0  +  (r4  -  r_4) 


Epc  -  Ec{ 0+) 
Jc0  =  -v  z  FCc0- 


D~ 


dC~ 


Avc  d  R 

The  potential  difference  Ec  —  £)(0+)  can  be  written  as 

£c-£c(0+)  =  ~^/hzc, 

^£hzc 


where 


(20) 

(21) 

(22) 

(23) 

(24) 


°hzc  —  rCact0CactFzMCO--  +  Tc6cFzMq-  >  (25) 

and  the  cathodic  external  current  as 


FSr, 


=  f 3  -  f-3  +  f4  -  T -4- 


(26) 


Notice  that  the  proposed  equations  have  a  correct  stationary  limit 
with  no  current:  in  fact,  few  mathematical  manipulations  lead 
to  the  classical  thermodynamical  Nernst  equation,  both  at  the 
anode  and  at  the  cathode  of  course: 

Ea-E„(  0+)  K-zK-2b  Cco^aoCH^OaO 

2F  —  In  __ - h  In  — zzrz: — - -  (27) 


RT 


K2K2b 


('aC)('H2aO 


Ec-Ec{ 0+)of  =  ^  F3F4  ^  (Cq2Co)(1/2)Cco2co 
ft  ~  K-3K-4  n  c;0 


(28) 


2.2.  The  space  charge  regions 


In  these  zones,  mass  conservation  equations  for  charged 
species,  both  negative  and  positive,  are  considered: 


dc±  d  j±  _ 

dt  +  dx 


(29) 


where,  again,  according  to  the  Planck-Nernst  formulation,  the 
flux  writes 


jf  —  —v±z±FC± 


dR 


D1 


zdC± 

Jw 


(30) 
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Notice  that,  again,  superscript  (“)  refers  to  the  carbonate  ion, 
while  superscript  (+)  refers  to  K+  and  Li+,  here  considered 
together  as  forming  a  unique  “equivalent”  cation. 

2.3.  The  axial  model  and  closing  equations 


In  that  part  of  the  electrolyte  stretching  from  one  electrode 
to  the  other,  the  mass  conservation  equation  for  a  hypothesised 
average  negative  charge  concentration  (in  fact,  electroneutrality 
writes  2C_ax  =  C+ ax),  namely 


9C~ax 
3 1 


+ 


94 


=  0 


(31) 


is  integrated  (from  ;tax  =  0  to  xax  =  /ax).  Substituting  the  flux 
using  the  Nernst-Planck  expression  yields  this  state  equation  in 
the  state  variable  C“ax: 


-  EC(8C)  -  Ea(8a )  C“ax  -  Cs 

C“ax  =  —2v~z~FC~ax  — — al  _  2D~  6c 


& 


Z2 

*ax 


2D 


-  Csa  C  ax 

llx 


(32) 


Finally,  the  circuit  is  closed  considering  the  Ohm  law  inside  the 
electrolyte  (this  is  part  of  the  axial  region  indeed),  i.e. 


*ele  —  ^ele" 


EC(8C)  ~  Ea(Sa) 

lax 


FSn 


the  external  electronic  current 


(33) 


text  —  (Ec  Ea)Gex  t  —  ic  —  ia 


(34) 


and  the  Kirchhoff  equation  for  potentials: 

Ec-Ea  =  (Ec  -  Ec( 0+))  +  (Ec( 0+)  -  EC(8C))  +  (EC(8C) 

-  Ea(Sa ))  +  (Ea(Sa)  -  £fl(0+))  +  (Ea(0+)  -  Ea).  (35) 


3.  Simulation  results 


This  section  reports,  first,  some  observations  about  parame¬ 
ter  estimation,  like  that  of  the  device  internal  resistance  and  of 
the  thickness  of  the  space  charge  regions,  and  then  some  sim¬ 
ulation  results  obtained  with  the  implemented  model,  focusing, 
in  particular,  on  the  stationary  power  versus  current  profile  and 
on  some  dynamical  characteristics.  For  simplicity,  here,  only 
an  elementary  MCFC,  instead  of  a  stack,  is  taken  into  account 
and  only  one  instance  (associated  to  the  whole  cell)  for  each 
electrochemical  variable  is  employed.  In  other  words,  the  ele¬ 
mentary  cell  is  considered  as  formed  by  a  unique  electrical  strip 
(corresponding  to  the  one  at  the  anode  and  cathode  inlet). 

Cell  internal  resistance  has  been  measured,  at  CESI,  on  a 
laboratory  stack,  through  a  short  circuit  experiment  [12]:  the  av¬ 
erage  result,  over  those  14  cells  with  700  cm2  useful  area  each 
and  over  different  hot  time  values,  was  0.7  m £2 .  The  value  ob¬ 
tained  with  a  1  |xs  short  circuit  simulation  is  0.698  mQ  (this  takes 
in  the  resistance  of  the  anodic  and  cathodic  metal  plates  and  that 
of  the  electrolyte,  but  not  that  of  the  bipolar  plate  and  of  the 
gas  distributor,  which  has  not  been  included  in  the  model  yet, 
but  which  is  presumably  negligible).  This  agreement  with  the 


experimental  value  contributes  to  support  model  reliability  (in 
conjunction  with  sound  parameter  calibration). 

A  simple  estimate  of  8a  and  of  8C  can  be  obtained  by  linearly 
approximating  the  ion  concentration  Galerkin  parabolas,  built 
on  the  three  concentrations  C~0,  C~{  and  C7  and  on  C~0,  C~t 

zi  2  c  ^ 

and  Cfc  respectively,  at  points  aO  and  a  4  and  at  cO  and  c\  re¬ 
spectively,  and  by  determining  the  vanishing  points  of  the  two 
straight  lines,  all  this  with  the  passing  of  time:  from  this  proce¬ 
dure,  it  appears  that  the  approximation  of  assuming  constant  8 , 
s  =  1(T8  m  is  rather  a  good  one. 

A  relevant  feature  in  the  cell  behaviour  is  the  dependence 
of  stationary  power  and  voltage  on  current:  model  reliability 
is  sustained  by  the  comparison  between  simulated  curves  such 
as  those  shown  at  the  top  of  Fig.  3  and  characteristic  curves 
found  in  the  literature  (compare,  for  example,  [13]).  The  fact 
that,  for  increasing  current  values,  power  starts  to  decrease  can  be 
explained  partly  as  a  kinetic  limit,  partly  as  due  to  concentration 
polarization:  in  fact,  for  high  current  values,  at  the  beginning  the 
problem  is  that  there  are  not  enough  activated  sites  to  sustain 
the  electrochemical  reactions  (as  shown,  again,  in  Fig.  3),  but 
then,  when  they  manage  to  be  produced,  the  carbonate  ion  is 
not  fast  enough  in  moving,  by  diffusion  most  of  all,  from  the 
anode  side  to  the  cathode  side,  to  feed  reactions  again;  therefore, 
the  potential  differences  Ea( 0+)  —  Ea  and  Ec  —  Ec( 0+)  (Fig.  3) 
decrease,  causing  the  external  potential  difference  Ec  —  Ea,  and 
therefore  power,  to  decrease.  This  explanation,  anyway,  could  be 
checked  and/or  refined  by  modifying  the  model,  in  particular  by 
inserting  a  more  truthful  pore  distribution,  instead  of  the  uniform 
one  adopted  so  far. 

As  already  remarked,  the  model  should  be  useful  also  to  in¬ 
quire  cell  dynamics;  therefore,  the  variations  exhibited  by  the 
electrochemical  variables  because  of  rather  fast  load  variations 
are  reported:  Fig.  4  shows  the  case  when  Gext  decreases,  from  a 
stationary  value  equal  to  1500  S,  to  1485  S,  in  1  s  (the  load  starts 
to  change  at  time  t  —  25  s  here).  Fig.  5  illustrates  an  example 
of  a  large  load  variation:  it  is  assumed  that  Gext  increases  from 
0  S  (at  time  t  =  0  s)  to  1500  S  in  10  s  and  then  sticks  to  the  final 
value  which  it  has  attained. 

4.  Conclusions 

A  dynamical  electrochemical  model  has  been  described  to 
study  the  fast  phenomena  involved  in  the  behaviour  of  a  MCFC. 
Simulation  results  coming  from  its  implementation  show  that  it 
is  able  to  reproduce  some  experimental  facts  like  power  decrease 
for  high  current  values. 

We  remark  that  the  model  can  be  integrated  with  a  thermo- 
fluid-dynamical  model  [5,6]  of  a  MCFC,  so  as  to  obtain  an  over¬ 
all  dynamical  description,  valid  for  both  low  and  relatively  high 
frequencies,  i.e.  for  time  scales  both  on  the  order  of  an  hour  or 
tens  of  seconds  and  on  the  order  of  a  second.  This  procedure  is 
presented  in  [14]. 
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